Quantum phase-space simulations of fermions and bosons 
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We introduce a unified Gaussian quantum operator representation for fermions and bosons. The 
representation extends existing phase-space methods to Fermi systems as well as the important case 
of Fermi-Bose mixtures. It enables simulations of the dynamics and thermal equilibrium states of 
many-body quantum systems from first principles. As an example, we numerically calculate finite- 
temperature correlation functions for the Fermi Hubbard model, with no evidence of the Fermi sign 
problem. 



Calculating the quantum many-body physics of inter- 
acting Fermi systems is one of the great challenges in 
modern theoretical physics. In even the simplest cases, 
first-principles calculations are made difficult by the com- 
plexity of the fermionic wave- function, manifest notori- 
ously in the Fermi sign problem. In previous quantum 
Monte Carlo (QMC) techniques, the sign problem ap- 
pears as trajectories with negative weights, which con- 
tribute to a large sampling errorQ, together with large, 
computationally intensive determinants. 

Fermion complexity issues appear in physical problems 
at all energy scales, from high-energy lattice QCD to the 
emerging area of ultra-cold atomic physics. Recent pio- 
neering experiments in ultra-cold Fermi gases are capable 
of investigating fermion many-body physics in regimes 
of unprecedented experimental simplicity. This situa- 
tion implies a substantial opportunity to develop and test 
novel first-principles theoretical methods for the investi- 
gation of correlations and dynamical effects. 

Here we present a phase-space method for simulating 
many-body boson0, El and fermion^, 0] systems, based 
on a Gaussian operator expansion. 

The method allows the treatment of dynamical and 
static problems at finite temperature. The expansion in 
the fermionic case represents pairs of Fermi operators. 
Since the pairs obey commutation relations, there are no 
anti-commutators causing sign problems, and no large 
determinant calculations as in some previous approaches. 

The method is illustrated using the finite tempera- 
ture Hubbard model, which is a well-known theory in 
condensed matter physics and high T c superconductors. 
The cases chosen have an acute sign problem using con- 
ventional QMC. The results are directly applicable to 
feasible experiments on ultra-cold fermions in an optical 
lattice 0. 

Like path- integral QMC, phase-space methods sample 
the many-body quantum density operator p. But rather 
than expressing the density operator in a position rep- 
resentation, one expands it in terms of an overcomplete 
basis of operators: 



where P( A , t) is a positive probability distribution, A 
is an overcomplete operator basis for the class of den- 
sity matrices being considered, and d A is the integration 
measure for the generalized phase-space coordinate A . 
It is the overcompleteness of the basis which allows a 
positive representation of any physical density matrix in 
terms of Gaussian operators. 

We define the operator basis A = f2A + A_ to be the 
product of Gaussian forms of bosonic ( A + ) and fermionic 
( A_) creation and annihilation operators, over M± 
modes respectively. Here Q is an additional weighting 
factor. If a is a row vector of M± annihilation operators, 
and the corresponding column vector of creation oper- 
ators, their commutation relations are: [ajt,oJ]q: = 5kj ■ 
We use =p to indicate bosons (upper sign) or fermions 
(lower sign). For brevity, we restrict the present discus- 
sion to number-conserving systems. The most general 
Gaussian operator is then a generalized thermal density 
operator with a complex covariance: 



A±(n) = |I±n 



exp 



a {21} -[I in]" 1 a 



(2) 



where the additional {1} in the exponent is bracketed 
to indicated that it only appears in the fermionic case. 
Normal ordering is defined as usual for Bose and Fermi 
systems, for example, : OjCij : = ±aj Sj = ±n,j . The 
M± x M± matrix n corresponds to a generalized thermal 
covariance. 

Using these Gaussian operators in the density operator 
expansion of Eq. (1), one finds that operator expectation 
values become weighted moments of the distribution P, 
denoted as (..)p. Thus the first- and second-order num- 
ber correlations are 



(n) = J ilnP(~t,t)d~t = (n) P , 

ajcii) = (nnnjj) P ± (n lj n ji ) 1 



(3) 



p(t) = / P(X,t)A(X)dX 



(1) 



As an illustration of the use of the unified represen- 
tation, consider the canonical distribution of a Bose or 
Fermi field. The thermal state at temperature T = 



2 



1/feeT can be cast into an imaginary time integro- 
differential equation: 



dp 

~d7 



i 



P(X,t)[H ,A(A)]+dA 



(4) 



To solve this, we first use identities derived in 0,0 that 
describe the action of operators on the density operator 
as derivatives on elements of the Gaussian basis. After 
integrating Eq (@J by parts, we arrive at the following 
mappings: 

no -> nP l^- (I±n) T l nP 
{on J 

pn -> nP-(^-n T l (I ± n) P . (5) 
on 



The matrix derivative is here defined as (d/dn)ij — 
d/dtiij . In the free-field case of H = "Y^u^k^ki we 
arrive at a first-order Fokker-Planck equation for the 
distribution function P. This leads to deterministic 
equations for the mode occupations of form dnk/dr = 
—uikTik (1 ± nit) , which can be integrated to get the well- 
known Bose-Einstein (Fermi-Dirac) distribution: 



exp(w fc T) =p 1 



(6) 



For systems of interacting particles, the unified rep- 
resentation gives nonlinear, stochastic phase-space equa- 
tions, which must be solved numerically. Due to the non- 
uniqueness of the expansion, a careful choice of identities 
is mandatory^ to keep the nonlinear equations stable. 
We term this choice a stochastic gauge in analogy to the 
gauge fields in QED, since it results in freely chosen fields 
in the resulting stochastic equations. The choice in the 
fermionic case is especially large, since the fermionic anti- 
commutation relations result in a non-unique algebraic 
form of the Hamiltonian. 

As an example, consider the Hubbard model, which 
is the simplest nontrivial model for strongly interacting 
electrons. It is an important system in condensed matter 
physics, with relevance to the theory of high-temperature 
superconductors The Hamiltonian is: 

PT(ni,n_i) = - (Uj + fJ-Sjj) n ijt<7 

-i^iE : (%i- s ^-i) 2: / 2 ( ? ) 



where n iha = aj^aj^ = {ri a }^. The coupling ty = t 
if the i,j correspond to nearest neighbour sites and is 
otherwise 0. The index a denotes spin (±1) and the 
indices i,j label lattice location, and s = U/\U\ = ±1. 
Traditional QMC methods for this problem have large 
sign problems in the repulsive case (s = 1)[H,E3- 



Because of the way we have written the interaction 
term in Eq. (6) (which constitutes a kind of fermionic 
gauge choice) , the Hubbard model maps to a set of real 
Stratonovich stochastic equations: 



^ = l{(I-n,)T(\ + IV Ti 2 )(I-n„)} 



(8) 



Here we have introduced the stochastic propagation ma- 
trix: 



T. 



(r) 



\ U \ S H i sn 33-° - n ii° + 7i } ■ ( 9 ) 



The real Gaussian noise Q r ' (r) is defined by the correla- 
tions 

(^VHjrV)) =2\U\S(r-r')S jjl S rrl . (10) 

The weights for each trajectory evolve as physically 
expected for energy- weighted averages, with dfl/dr = 
— OP^(ni, n_i). Because the equations for the phase- 
space variables n, Ji( , are all real, the weights of all tra- 
jectories will remain positive. Thus the traditional man- 
ifestation of the sign problem is avoided, as there is no 
'deterioration of the sign' from averaging over positive 
and negative weights. Furthermore, the mapping to real 
phase-space produces a stable set of equations, and thus 
there is no need to invoke additional gauge choices. 

There is, however, the issue of spreading weights, which 
becomes serious for large lattice sizes and long simula- 
tions times. Physical quantities are weighted averages: 



N,-, 



(A(r)) p = V {p) (r)# (r)/ £ fi« ( T ) 



(11) 



0=1 



0=1 



where N p is the total number of paths in the sample. 
A large spread in the weights makes a straightforward 
average very inefficient, as most paths in the ensemble 
may end up contributing very little to the final result. 
To increase efficiency, we instead use a simple branching 
algorithm adapted from Green's function Monte Carlo 
methods 0, in which low- weights paths are deleted and 
high-weight paths are cloned, according to the rate: 



= Integer £ + n (p) /n 



(12) 



where £ is a random variable uniformly distributed on 
[0, 1] and where f2 is a reference weight, which is adapted 
to keep the number of paths N p under control. At 
branching, the weights are equalised and thereafter the 
clones evolve independently with spreading weights. To 
avoid biasing, the branching must occur sufficiently often 
to limit the diversity of weights at the branching times. 
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For the results presented here, the branching algorithm 
is sufficient to control sampling error - other situations 
may require the use of more sophisticated importance 
sampling methods 0- 

The stochastic phase-space equations are simulated by 
a robust semi-implicit algorithm [l^. with an adaptive 
stepsize to overcome stiffness. Unlike Projector QMG 
methods, the Gaussian phase-space method can calcu- 
late any correlation function, at any temperature. Un- 
like Path Integral QMC, a single run generates results 
for a range of temperatures: longer simulation times cor- 
respond to lower temperatures. Strictly speaking, zero- 
temperature results are obtained only in the limit of long 
simulation times. In practice, however, one only has to 
run the simulation until the relevant correlation functions 
have plateaued. 

Precision is of course limited by sampling error, but 
this can be reduced by several means. For example, 
one can a) include more trajectories in the sample, b) 
employ a more sophisticated branching/importance sam- 
pling technique to reduce the spread in weights, and c) 
make a better 'stochastic gauge' choice to obtain phase- 
space equations with smaller sampling error for the cor- 
relations to be calculated. 

Typical results for a 16 x 16 lattice are shown in Figs 
(1) and (2), which plot the energy E and second-order 
correlation function g2, respectively, for different chemi- 
cal potentials. The estimation of sampling error shown 
in the figures assumes independent samples, for simplic- 
ity. While this is liable to underestimate the error, espe- 
cially for <?2 where there is also spatial averaging, it does 
indicate the approximate dependence on temperature. 
In particular, the sampling error remains well-controlled 
throughout the simulation, even away from half filling, 
where there is known to be a sign problem. A more de- 
tailed sampling error analysis will be given elsewhere. 

In conclusion, we have presented a unified operator 
representation that is able to represent arbitrary physi- 
cal states of bosons and fermions. By use of this repre- 
sentation, non-interacting systems can be mapped to de- 
terministic phase-space equations, whereas systems with 
two-body interactions can be simulated by use of stochas- 
tic sampling methods, provided a suitable gauge is chosen 
to eliminate any boundary terms. For the example of the 
Hubbard model, we show how the thermal equilibrium 
problem can be mapped to a set of real, stable phase- 
space equations with positive weights. Bosonic problems 
can be solved in similar manner, and the method can also 
be used to simulate dynamics (although typically with an 
increased sampling error). Thus the one, unified method 
can solve both fermionic and bosonic problems, which 
makes it well suited to simulating Bose-Fermi mixtures, 
and to studying the BEC/BCS crossover. 
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Figure 1 : Total energy E per site versus inverse temperature r 
for a 16 x 16 2D lattice for chemical potentials /i = 2 (solid), 
fi — 1 (dashed) and fi — (dot-dashed). Curves without 
crosses give the number of particles per site for fi = 1 (dashed) 
and \x = (dot-dashed). U = 4, t = 1, and 100 paths initially. 
Dotted curves give an approximate sampling error. 
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Figure 2: Second-order correlation function gi = 
) / "m.i ) versus inverse 

temperature r for a 16 x 16 2D lattice for chemical poten- 
tials jj, — 2 (solid), fj, — 1 (dashed) and /j, = (dot-dashed). 
U = 4, t = 1, and 100 paths initially. Dotted curves give an 
approximate sampling error. 



[1] D. M. Ceperley, Rev. Mod. Phys. 71 (1999) 438. 
[2] J. F. Corney, P. D. Drummond, Phys. Rev. A 68 (2003) 
063822. 

[3] P. D. Drummond, P. Deuar, and K. V. Kheruntsyan, 

Phys. Rev. Lett. 92, (2004) 040405. 
[4] J. F. Corney, P. D. Drummond, Phys. Rev. Lett. 93 

(2004) 260401. 
[5] J. F. Corney, P. D. Drummond, cond-mat/0411712. 
[6] L. Pezze, L. Pitaevskii L, A. Smerzi, S. Stringari, G. Mod- 

ugno, E. de Mirandes, F. Ferlaino, H. Ott, G. Roati, 

M. Inguscio, Phys. Rev. Lett. 93, 120401 (2004). 



4 



[7] P. Douar, P. D. Drummond, Phys. Rev. A 66 (2002) 
033812. 

[8] W. von der Linden, Phys. Rep. 220 (1992) 53. 
[9] R. R. dos Santos, Brazilian Journal ol Physics 33 (2003) 
36. 

[10] W. Fettes, I. Morgenstern, Computer Physics Communi- 



cations 124 (2000) 148. 
[11] N. Trivedi, D. M. Ceperley, Physical Review B 41 (1990) 
4552. 

[12] P. D. Drummond, I. K. Mortimer, J. Computat. Phys. 
93 (1991) 144. 



